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In recent years, a second fluid-fluid phase transition has been reported in several materials at 
pressures far above the usual liquid-gas phase transition. In this paper, we introduce a new model 
of this behavior based on the Lennard- Jones interaction with a modification to mimic the different 
kinds of short-range orientational order in complex materials. We have done Monte Carlo studies of 
this model that clearly demonstrate the existence of a second first-order fluid-fluid phase transition 
between high- and low-density liquid phases. 



PACS number: 64.70.Ja 

The most common example of a first-order phase tran- 
sition is that between a liquid and a gas, such as boil- 
ing water. On the other hand, while many transitions 
between different solid phases of homogeneous materials 
are also well known, it is only relatively recently that ev- 
idence of a second fluid-fluid phase transition has been 
found. In fact, liquid-liquid phase transitions (LLPT) 
have been suggested in liquid S, Ga, Se, Te, I2, Cs, and 

Stell and Hemmer [|2|,[5| showed the existence of LLPT 
in a one dimensional model with a softened hard core po- 
tential and a long-range negative attraction. Their work 
was later studied in more detail by Franzese, Malescio, 
et. al. Q and Sadr-Lahijany, Scala, et. al. H. 

Mitus, Patashinskii and Shumilo Q proposed LLPT 
in molten salt at high pressure based on a phenomeno- 
logical model. Ferraz and March suggested a similar 
LLPT in carbon, with indirect experimental evidence be- 
ing found by Togaya 0. Glosli and Ree || published 
results of a first-order liquid-liquid phase transition in 
molten Carbon between two thermodynamically stable 
liquid phases. Extensive computer simulations on models 
of water have supported the existence of a LLPT in the 
metastable region [p[-[f6[. Experimental results support- 
ing the evidence of liquid-liquid phase transitions in wa- 
ter have also been found fl7H19|. Katayama and Mizu- 
tani et. al. p0| found a liquid- liquid phase transition in 
molten phosphorus using x-ray diffraction. 

Our objective is to investigate the general phenomenon 
instead of studying LLPT for any particular substance. 
We have developed a simple model that exhibits a tran- 
sition between high- and low-density liquids at high pres- 
sure. 

The behavior of our model is constructed to be similar 
to behavior seen in simulations of real substances such as 
water and carbon but without introducing the complexity 
of having to simulate molecular orientations as in water 
and carbon. 

To mimic the effects of local ordering, we have rep- 
resented the different relative local orientations of the 
molecules with a spin-one-half variable. The interactions 
between particles with the same spin are given by the 
original Lennard- Jones expression, while the interactions 
between particles with opposite spin are purely repulsive. 
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Note that there are different values of a for like and unlike 
spins. This is an important feature of the model, and 
some properties, including the symmetries of the solid 
phases, are sensitive to the relative values of ai and a u . 
The LLPT occurs when a u is smaller than 01 , so that by 
re-orienting the spins, the particles are capable of forming 
different co-ordination numbers and local structures. We 
have performed most of our investigations for the case in 
which the ratio is 1/2. 

Our model can be modified to couple to a fictitious ex- 
ternal magnetic field. If we define the magnetic moment 
as the sum of all spins, then the Hamiltonian is given by, 
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h is the fictitious magnetic field which is set to zero in 
our simulations and er^ = ±1 is the direction of the spins 
of individual particles. 

We used the units T* = k B T/e, P* = afk B P/e. We 
have set the potential cutoffs at 3er; . 

We performed Monte Carlo simulations of the two- 
dimensional version of our model in various ensembles. 
Although our simulations are done in two dimensions, 
our model is not limited to two dimensions. The simula- 
tions turned out to be rather difficult, and it was neces- 
sary to extend the usual techniques to improve efficiency. 
However, we did find clear evidence of an LLPT at high 
pressures. We are able to map out both PT and Pp phase 
diagrams. 

For our simulations of fluids, we confined the particles 
to a square with periodic boundary conditions. For those 
simulations that were extended to include solid phases, 
we used parallelograms to allow us to vary the angle of 
the boundary conditions, as well as the volume of the 
container. The optimization of the Monte Carlo simula- 
tion step sizes for individual particle motion, changes of 
the volume of the box, and changes of the angle of the 
parallelogram were dynamically optimized using the Ac- 
ceptance Ratio Method |^l[ . Metropolis flips to maintain 
equilibrium for the spins associated with each particle 
were also carried out. 
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Volume changes on the dense fluid turned out to be 
rather inefficient because <f>-\i increases very rapidly at 
short distances. We solved this problem by introducing 
a cluster Monte Carlo move. The clusters are formed by 
creating bonds between particles with probability, 
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and then changing the size of the box by rescaling the 
locations of the centers of mass of the clusters. This 
cluster move is extremely effective, with improvements 
in the acceptance ratio by factors of up to 10 10 . 

The probability of accepting a proposed cluster move 
is given as follows. Let £ be the set of all edges which 
joins any two particles which are less than 3ct; apart, and 
let B be the subset of £ consisting of bonds formed with 
the probability P(r), where r is the edge length joining 
the two particles. The probability of forming B is, 
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Hence the acceptance probability is, 
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Another new move that has proven extremely effective 
is to form clusters of nearby particles that are less than 
1/ \/ psin(n/3) apart and attempt to flip the spins of all 
particles in a cluster. 

For canonical ensemble simulations using two boxes 
in equilibrium with each other, the total volume is con- 
served pq] . An additional Monte Carlo move is intro- 
duced to transfer particles between the boxes. 

Thermodynamic quantities were calculated from simu- 
lations in the constant pressure ensemble over relatively 
large pressure and temperature domains using the Mul- 
tiple Histogram Method |p2| . The coexistence curve is 
mapped out by tracing the ridge line of the isothermal 
compressibility. The position of the tricritical point was 



estimated using finite size scaling. The Pp diagram was 
determined by combining other data with the results of 
simulations using the canonical ensemble (23| . 

The phase diagram for our model is shown in figure 
(|l|) . At low pressures its behavior is virtually identical to 
that of the usual Lennard- Jones model. At higher pres- 
sures, the coexistence curve for the LLPT is shown pro- 
jecting into the fluid region of the phase diagram. This 
line of first-order phase transitions actually ends in a tri- 
critical point, which separates it from a line of critical 
points associated with the ordering of the "spins" used 
to define the model. This line joins the liquid-gas and 
liquid-liquid first-order coexistence curves at the tricrit- 
ical points. However, it is only our access to the "mag- 
netic" degrees of freedom in the model that make it ob- 
vious that the LLPT curves end in tricritical rather than 
critical points. If we did not have access to the magnetic 
degrees of freedom, rather careful measurement of the ex- 
ponents characterizing the divergences at the end of the 
line of first-order transitions would be necessary to distin- 
guish the two cases. It is not completely clear whether a 
fictitious ordering field in a more elaborate model would 
not have the same consequence of producing tricritical 
points. 

The apparent discontinuity of the coexistence curves 
between the high-density solid, the low-density solid and 
the low-density liquid is an artifact of the accuracy of our 
simulations. These lines are obtained from the Multi- 
ple Histogram extrapolations in opposite directions. The 
fact that these two lines do not meet exactly, but termi- 
nated very close to each other shows the accuracy of our 
simulations. 

To determine the liquid-liquid coexistence curve in the 
PT diagram, simulations were done with 160 particles 
at constant pressure, using larger systems to check the 
results. Sets of histograms were collected at T*=0.53, 
0.55, and 0.565 over a range of pressures covering across 
the coexistence region. The system was equilibrated for 
30,000 MCS/P before beginning to take data for another 
300,000 MCS/P. 

To determine the tricritical point for the liquid-liquid 
phase transition, additional simulations for 160, 240, 320, 
and 480 particles were done at constant pressures at 
T*=0.55 and 0.565. The Monte Carlo simulations used a 
total of 300,000 MCS/P for the 160 particle systems and 
typical equilibration runs of 30,000 MCS/P. For larger 
systems, longer runs were used. For the 480 particle sys- 
tem, a total of 1,200,000 MCS/P were used with typical 
equilibration runs of 100,000 MCS/P. 

To determine the liquid-solid coexistence curve, sim- 
ulations were performed at P*=0.5, 0.6, 0.78 and 0.9 
over a range of temperatures covering the solid-liquid 
coexistence region. The system was equilibrated for 
60,000 steps per particle before collecting data for an- 
other 600,000 steps per particle. 

The phase transition between high- and low-density 
phases is strongly first-order. Figure (^a,b) shows snap 
shots of high- and low-density liquids near the coexis- 
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tence curve for a system with 240 particles. The low- 
density liquid is dominated by parallel nearest-neighbor 
interactions with a stronger core repulsion and with co- 
ordination number six. The high-density liquid is dom- 
inated by an anti-parallel nearest-neighbor interaction 
that allows the particles to come closer together and with 
co-ordination number three. The differences in local or- 
derings for low- and high-density liquids stabilize the liq- 
uids and create the possibility of liquid-liquid phase tran- 
sitions. 
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The volume-energy probability density function (fig- 
ure (Q)) has a saddle point that is typical of a first-order 
transition. The coexistence curve for the LLPT has a 
negative slope, as shown in the PT diagram in figure (Q), 
reflecting the higher entropy of the high-density liquid. 
The sharp bend of the liquid-solid coexistence curve at 
the triple point found from the multiple-histogram anal- 
ysis is consistent with our observation of the liquid-liquid 
phase transition. 

To find the tricritical point, we use the fact that peak 
values of the isothermal compressibility grow as 0(N) 
at the coexistence curve and as £>(iV a925 ) pi|-|29|] at the 
tricritical point, assuming that the transition is in the 
expected two-dimensional Ising class. Simulations were 
done with N=160, 240, 320 and 480 near the tricriti- 
cal point. Peak values of the isothermal compressibility 
are then plotted against system sizes for various tem- 
peratures. Figure (pj) shows the size dependence of the 
isothermal compressibility. Although it is hard to deter- 
mine the exact location of the tricritical point, the data 
suggests that it is located at T* = 0.56 ± 0.01. 
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FIG. 1. (a) PT diagram. Solid lines show coexistence 
curves obtained from our simulations with 160 particles, o 
marks the calculated solid-solid transition point at T=0. Dot- 
ted lines show our estimate for the coexistence curves. The 
thick dashed line represents the locus of Curie points for the 
second-order magnetic transitions. |m| and |m*| are the mag- 
netic and anti-magnetic order parameters respectively, (b) 
PT diagram in the region of interest, the circle marks the 
tricritical point and dotted line shows the peak of the isother- 
mal compressibility beyond tricritical point. Solid lines show 
coexistence curves obtained from our simulations with 160 
particles. The dashed line shows the extent of curve shifting 
due to the finite size effect. The upper dashed line corre- 
sponds to a system with 480 particles and the lower dashed 
line corresponds to a system with 320 particles. 



FIG. 2. Magnetic susceptibility plotted against T* 
P* = 0.5 and h = 0. 
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Figure g shows the growing divergence of the mag- 
netic susceptibility for 160, 240 and 320 particles. These 
magnetic susceptibility are plotted at P* = 0.5 which is 
far away from the liquid-liquid coexistence curve and the 
gas-liquid coexistence curve. 

The p-T diagram (figure (||)) was mapped out using 
simulations in the canonical ensemble using dual boxes 
to simulate coexistence without the inconvenience of an 
interface between the phases [^3| . The low-density liquid- 
gas coexistence region is only slightly lower than that of 
pure Lennard-Jones ]3(|[3l| due to the small effect of the 
differences in the models at low densities. 
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(a) 
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FIG. 3. (a,b) Snapshots of a system of 240 particles in the 
low-density(a) and high-density(b) liquid state coexisting at 
T*=0.55, P*=0.75, p*=0.833 and 1.25 respectively. 





FIG. 4. Probability density function for 240 particles coex- 
isting at T*=0.54 and P*=0.75. 
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FIG. 5. p-T diagram. Points obtained by simulation 
using the canonical ensemble with 160 particles in two 
boxes. Dashed lines show the coexistence region of the pure 
Lennard- Jones system reported by Barker, Henderson, Abra- 
ham and Phillips, Bruch, Murphy. Dotted lines show the 
coexistence region of our model. All lines are drawn to guide 
the eye. 
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FIG. 6. Plot of peak values of isothermal compressibility 
vs system size. The tricritical point is near T*=0.56. 

Figure(f?]a) shows a snapshot of the high-density crys- 
talline state, which has three-fold symmetry and zero 
magnetic moment. The low-density crystalline state (fig- 
ure(0b)) is hexagonal close packed with uniform spin. 
At zero temperature, we can calculate the location of 
the transition between the two solid phases to arbi- 
trary precision. Its location is determined to be P* — 
0.534819±1CT 7 , with p* = 1.82646 and p* = 0.943122 for 
the high- and low-density solids respectively. For com- 
parison, p* = 0.934721 at T* = and P* = 0, so the 
low-density solid phase changes its density very little up 
to the boundary of the high-density phase. 
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(b) 

FIG. 7. (a,b) The high- density (a) and low-density(b) solid 
states, notice a vacancy in the low-density solid state. 

It was found that symmetry of the solid phase is highly 
dependent on the ratio t = a u jai. Figure (|J) shows a 
solid with four-fold symmetry obtained from a simulation 
with r = 0.6. We believe that the full range of solid 
phases is quite rich for this model. 




FIG. 8. Four fold solids form using a model with r = 0.6. 

We have developed a relatively simple model that 
demonstrates a liquid-liquid phase transition between 
high- and low-density phases. Although our simulations 



are two dimensional, our model is not restricted to two 
dimensions and we see no reasion to believe that a three 
dimensional simulation with our model would produce 
significantly different behavior. By comparing the be- 
havior of our model with the properties of real systems, 
we hope to learn which properties are generic and which 
depend on details of a particular material. One immedi- 
ate point of interest is the negative slope of the liquid- 
liquid coexistence line in our model, which reflects the 
high entropy in the high-density phase. This feature is, 
indeed, found in most materials that exhibit an LLPT. 
However, this does not appear to be universal, since the 
coexistence curve of molten carbon reported by Glosli 
and Ree || nas a positive slope. At the present time, 
this difference in materials properties is not understood. 
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